% fig12_a1_adaptive.m ——— 自适应网络下 Fig.12(a1)：\tau=40 时的平均膜电位 y
clear; clc; close all;

N_vals = 1:50;
k_vals = linspace(0.1,5,50);
tau    = 40;

Nn = numel(N_vals);
Kn = numel(k_vals);

Y_mean = nan(Nn, Kn);

parfor iN = 1:Nn
    N = N_vals(iN);
    row = nan(1, Kn);
    for ik = 1:Kn
        k0 = k_vals(ik);
        y = simulate_node_adaptive(N, k0, tau);
        row(ik) = mean(y);
    end
    Y_mean(iN,:) = row;
end

delete(gcp('nocreate'));

figure;
imagesc(k_vals, N_vals, Y_mean);
axis xy;
set(gca,'XTick',linspace(k_vals(1),k_vals(end),6));
xlabel('k','FontSize',12);
ylabel('N','FontSize',12);
title('Fig.12(a1) — \langle y\rangle at \tau=40 (adaptive)','FontSize',14);
h=colorbar;
ylabel(h,'\langle y\rangle','FontSize',12);


function y_vec = simulate_node_adaptive(N, k0, tau)
% simulate_node_adaptive  对 N×N 网络在时间 τ 处做一次快照
%   y_vec = simulate_node_adaptive(N, k0, tau)
%   返回展平后的膜电位向量 y_vec（长度 N^2）
%
% 采用 Euler 法积分，并按论文式 (13),(15) 更新耦合强度 k_ij 和参数 alpha_ij。

    % —— 全局模型参数 —— 
    a     = 0.5;    b    = 0.215;
    beta  = 0.2;    A    = 0.95;
    omega = 0.4;    mu_s = 0;

    % —— 自适应增长律参数 —— 
    r      = 0.0007;    % 耦合强度增益
    sigma2 = 0.005;     % alpha 增益
    eps1   = 1e-5;      % alpha 阈值
    eps2   = 1e-5;      % k 阈值

    % —— 积分设置 —— 
    dt   = 0.01;
    T0   = 200;                 % 丢弃瞬态
    Nt0  = round(T0/dt);
    Nt   = round(tau/dt);

    % —— 初始条件 —— 
    X = zeros(3, N, N);
    X(1,:,:) = 0.01;
    X(2,:,:) = 0.02;
    X(3,:,:) = 0.02;

    K     = k0 * ones(N,N);
    Alpha = 1.8 * ones(N,N);
    Alpha(15:35,15:35) = 1.795;

    H_prev = zeros(N,N);

    %—— 辅助：耦合计算 ——
    function W = coupling(Y)
        Yp = padarray(Y, [1 1], 'replicate');
        W  = zeros(N,N);
        for ii=1:N
            for jj=1:N
                y0 = Yp(ii+1,jj+1);
                nb = Yp(ii  ,jj+1) + Yp(ii+2,jj+1) + ...
                     Yp(ii+1,jj  ) + Yp(ii+1,jj+2);
                W(ii,jj) = K(ii,jj)*(nb - 4*y0);
            end
        end
    end

    %—— 积分循环 —— 
    for step = 1:Nt
        Ymat = squeeze(X(2,:,:));
        W    = coupling(Ymat);

        for i=1:N
            for j=1:N
                x = X(1,i,j); y = X(2,i,j); z = X(3,i,j);
                k_ij   = K(i,j);
                alpha_ij = Alpha(i,j);

                dx = -alpha_ij*((y-x)-(y-x)^3/3 + (a+3*b*z^2)*x) + W(i,j);
                dy =      (y-x)-(y-x)^3/3 - beta*y + mu_s ...
                      + A*cos(omega*step*dt) + W(i,j);
                dz =      x + W(i,j);

                X(1,i,j) = x + dt*dx;
                X(2,i,j) = y + dt*dy;
                X(3,i,j) = z + dt*dz;
            end
        end

        if step > Nt0
            % 计算当前能量
            X_flat = reshape(X, 3, []);
            H_vec  = compute_energy(X_flat, mean(Alpha(:)), a, b);
            H_cur  = reshape(H_vec, N, N);

            for i=1:N
                for j=1:N
                    % 更新 Alpha
                    if H_cur(i,j) - H_prev(i,j) > eps1
                        Alpha(i,j) = Alpha(i,j) + sigma2*Alpha(i,j)*dt;
                    end
                    % 更新 K
                    neigh = H_cur(max(i-1,1):min(i+1,N), max(j-1,1):min(j+1,N));
                    if H_cur(i,j) - mean(neigh(:)) > eps2
                        K(i,j) = K(i,j) + r*dt;
                    end
                end
            end
            H_prev = H_cur;
        end
    end

    Y_final = squeeze(X(2,:,:));
    y_vec   = Y_final(:);
end

